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We perform a density-matrix renormalization group (DMRG) study of the S = | Heisenberg an- 
tiferromagnet on the kagome lattice to identify the conjectured spin liquid ground state. Exploiting 
SU(2) spin symmetry, which allows us to keep up to 16 000 DMRG states, we consider cylinders 
with circumferences up to 17 lattice spacings and find a spin liquid ground state with an estimated 
per site energy of —0.4386(5), a spin gap of 0.13(1), very short-range decay in spin, dimer and 
chiral correlation functions and finite topological entanglement 7 consistent with 7 = log 2 2, ruling 
out gapless, chiral or non-topological spin liquids in favor of a topological spin liquid of quantum 
dimension 2, with strong evidence for a gapped topological Z2 spin liquid. 
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A pervasive feature of physics is the presence of sym- 
metries and their breaking at low energies and tempera- 
tures. It would be an unusual system in which at T = 
(quantum) fluctuations are so strong that all symme- 
tries remain unbroken in the ground state. In mag- 
netic systems such a state is dubbed a quantum spin 
liquid (QSL) [I] and is most likely to occur if fluctu- 
ations are maximized by low-dimension, low-spin and 
strong geometrical frustration; the search for a QSL has 
thus focused on frustrated S = \ quantum magnets in 
two dimensions. The Heisenberg antiferromagnet on the 
kagome lattice[2] (KAFM) is a key candidate, described 
by the S = \ model 

« = £3-4 a) 

with (i, j) nearest neighbors. 

Experimentally, the focus is on the Herbertsmithite 
ZnCu 3 (OH) 6 Cl 2 , modeled by Q on a kagome lattice 
with additional Dzyaloshinskii-Moriya interactions [3]- It 
is thought that the ground state is a spin liquid[4*HTU]. 
with no onsite magnetization [HJ [TT] and no spin gap [111 - 
IT4] within very tight experimental bounds. 

On the theoretical side, the kagome model of Eq. ([!]) re- 
mains a formidable challenge. While all proposed ground 
states show no magnetic ordering, they can be classified 
by whether they break translational invariance or not. 
The former type of ground state, a valence bond crys- 
tal (VBC), was pioneered by Marston|15j. The emerg- 
ing proposal was that of a "honeycomb VBC" (HVBC) 
with a hexagonal unit cell of 36 spins [TSIED] sharing in 
dimer-covered hexagons and a sixfold "pin wheel" at the 
center. On the other hand, a multitude of QSL states 
were proposed [2TH32] . Proposals for a QSL ground state 
include a chiral topological spin liquid (2TJ [22l [331 131], a 
gapless spin liquid 23-26], and various Z 2 spin liquids[2"?l- 



[30] with topological ground state degeneracy. 

In the past, numerical methods failed to resolve the 
issue conclusively. Quantum Monte Carlo faces the sign 
problem. Sizes accessible by exact diagonalization H |3"51 - 
|4"5] are currently limited to 48 sites. Other approaches 
diagonalized the valence bond basis or applied the con- 
tractor renormalization group (CORE) method, or the 
coupled cluster method (CCM) [4"M5"5] . The multi-scale 
entanglement renormalization ansatz (MERA) [56 found 
the VBC state lower in energy than the QSL state re- 
ported in an earlier density-matrix renormalization group 
(DMRG) study of tori up to 120 sites [31"] . 

Recently, strong evidence for a QSL was found in a 
large-scale DMRG study [57] considering long cylinders 
of circumference up to 12 lattice spacings. Ground state 
energies were substantially lower than those of the VBC 
state and an upper energy bound substantially below the 
VBC state energy was found; the ground state, having 
the hallmarks of a QSL, was not susceptible to attempts 
to enforce a VBC state. As to the type of QSL, [57] 
did not provide direct evidence for a Z 2 topological QSL. 
This has sparked a series of papers trying to identify the 
QSL [23 [21 [31 [SI EE], where again chiral spin liquids 
and gapless U(l) spin liquids were advocated and a clas- 
sification of Z 2 spin liquids achieved. At the moment, 
the issue is not conclusive. 

Here we study the KAFM using DMRG [5M)T] . in 
the spirit of [57 . DMRG is a variational method 
in the ansatz space spanned by matrix product states 
(MPS) which allows it to find the ground state of one- 
dimensional (ID) systems efficiently even for large sys- 
tem sizes. It can also be applied successfully to two- 
dimensional (2D) lattices by mapping the short-ranged 
2D Hamiltonian exactly to a long-ranged ID Hamilto- 
nian [57j [62l466j . DMRG cost scales roughly exponen- 
tially with entanglement entropy, such that area laws 
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TABLE I. Ground state energy per site (E/N) and gaps for 
L — oo cylinders (circumference c). Errors are from extrapo- 
lation; comparisons are with Ref. [57] except for the tori. 
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FIG. 1. Bulk energies per site. Lengths are in units of lattice 
spacings. The HVBC result |18l 119). and the upper bounds 
of MERA [56] and DMRG [57] apply directly to the thermo- 
dynamic limit; 2D estimates are extrapolations. 



limit system sizes, and DMRG favors open boundary con- 
ditions (OBCs) over preferable periodic boundary condi- 
tions (PBCs). The conventional compromise [57], taken 
also by us, is to consider cylinders, i.e. PBCs along the 
short direction (circumference c) and OBCs along the 
long direction (length L) where boundary effects are less 
important. Cost is dominated exponentially by circum- 
ference c. We use two different ID mappings (labeled as 
XC and YC plus cylinder size) [57] to check for undesired 
mapping dependencies of the DMRG results. Instead of 
earlier Abelian U(l) DMRG with up to 8 000 ansatz 
states, we employ non-Abelian SU(2) DMRG [SO EH] 
based on irreducible representations corresponding to 
16 000 ansatz states in a U(l) approach. This has crucial 
advantages: available results can be verified with much 
higher accuracy. The circumference of the cylinders can 
be increased by almost 50 % from 12 to 17.3 lattice sites 
(up to 726 sites in total) , strongly reducing finite size ef- 
fects; we also consider tori of up to 108 sites. We can 
eliminate the spin degeneracy that necessitates pinning 
fields in U(l)-symmetric simulations and avoid artificial 
constraints in gap calculations, making them more accu- 
rate and reliable. We also present results on spin, dimer 
and chiral correlation functions, the structure factor and 
topological entanglement entropy. All data agree with a 
gapped non-chiral Z 2 spin liquid; other QSL proposals 
for the KAFM are inconsistent with at least one of the 
numerical results. 

.Energies. -Energies for cylinders of fixed c and L are ex- 
trapolated in the truncation error of single-site DMRG 
[70] ; bulk energies per site are extracted by a subtrac- 
tion technique \§jg\ and extrapolated to L — > oo. Re- 
sults for various ID mappings and c are displayed in Ta- 
ble |TJ We also show the spin (triplet) gap to the 5 = 1 
spin sector. We confirm and extend earlier results [57] , 
At 16 000 states, DMRG is highly accurate; negligible 
changes in energy for substantially larger c support that 



the thermodynamic limit energy is found, which we place 
at -0.4386(5) (Fig. [TJ. Similar to Ref. [ST] we find the 
energy to be significantly below that of VBC states and 
no trace of a VBC in the correlation patterns. Except 
for the edges, bond energies are fully translationally in- 
variant. All results are consistent with strict variational 
upper bounds obtained without extrapolations from in- 
dependent DMRG calculations for infinitely long cylin- 
ders using the iDMRG variant [71], which are below the 
VBC energies. 

On the issue of a spin (triplet) gap [35J [72], Yan et 
al. [57) argue in favor of a small, but finite spin gap. 
SU(2) DMRG computes the 5 = 1 state directly and 
more efficiently; boundary excitations are excluded by 
examining local bond energies. We find the spin gap 
(Table |l] and Fig. |J to remain finite also for cylinders 
of large c. Whereas the results for small c agree with 
the 5 = 1 state energies and gaps reported in [57J, they 
display significant differences for larger c, perhaps due 
to the more complex earlier calculation scheme. SU(2)- 
invariant results evolve more smoothly with c, allowing 
a tentative extrapolation to a spin gap = 0.13(1) in 
the thermodynamic limit. Size dependence is small, in 
line with very short correlation lengths. The finite spin 
gap contradicts conjectures of a U(l) or other gapless 
spin liquids. For the calculation of the singlet gap found 
to be finite in Ref. [57], SU(2) DMRG does not offer a 
significant advantage to be reported here. 

Correlation Functions. -For all cylinders, we find an 
antiferromagnetic spin-spin correlation function (5; • Sj) 
along different lattice axes with almost no directional de- 
pendence. Exponential fits with a very short correlation 
length of £ ~ 1 (Fig. 4(a) ) were consistently better than 



power law fits, in agreement with a spin gap. This is not 
consistent with an algebraic spin liquid [23] , where the 
correlations are predicted to decay according to a power 
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We also consider the static spin structure factor S(q) 
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5 lq ' { - ri ~ Tj > (Si-Sj), gin units of basis vectors (61,62) 
of the reciprocal lattice. The spectral weight is concen- 
trated evenly around the edge of the extended Brillouin 
zone, with not very pronounced maxima on the corners of 
the hexagon. Results for large cylinders agree well with 
ED results for tori up to 36 sites [33]. All our S(q) are 
in accordance with the prediction for a Z 2 QSL|27j. 

We also find antiferromagnetically decaying, almost 
direction-independent dimer-dimer correlations for which 
again an exponential fit is favored (Fig |4(b)j ), in agree- 
ment with a singlet gap. Our data do not support the 
algebraic decay predicted [33] for an algebraic QSL. 

Chiral correlation functions [30] (CV,-fcC; mn ) = (Si ■ 
(Sj x Sk) ■ Si ■ (S m x S n )), where the loops considered 
are elementary triangles, did not show significant cor- 
relations for any distance or direction and decay expo- 
nentially (Fig. [5]), faster than the spin-spin correlations. 
Expectation values of single loop operators Cijk vanish, 
as expected for finite size lattices. Chiral correlators for 
other loop types and sizes decay even faster. Our findings 
do not support chiral spin liquid proposals [2TJ [531 133] • 

Topological Entanglement Entropy.-To obtain direct 
evidence regarding a topological state, we consider 
the topological entanglement entropy [73H75] . For the 
ground states of gapped, short-ranged Hamiltonians in 
2D, entanglement entropy scales as S ~ c, if we cut 
cylinders into two, with corrections in the case of topo- 
logical ground states [75]. We examine Renyi entropies 
S a = (1 — a) -1 log 2 trp a , < a < 00, where p is a subsys- 
tem density matrix. Scaling is expected as S a ~ rjc — 7 
where r\ is an a-dependent constant. 7, the topologi- 
cal entanglement entropy, is independent of a [77H79] 
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FIG. 2. Plot of the bulk triplet gap for infinitely long cylinders 
versus the inverse circumference c in units of inverse lattice 
spacings with an empirical linear fit to the largest cylinders, 
leading to a spin gap estimate of 0.13(1). 



(a) XC16 system (196 sites) 



(b) 27-site torus system 



FIG. 3. Two static structure factors S(q); k x , ky in units of 
reciprocal lattice basis vectors. Results are independent of 
the choice of ID mapping (not shown). 
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FIG. 4. Log-linear plots of the absolute value of the |4(a)| 
spin-spin and 4(b) dimer-dimer correlation functions versus 
the distance x~=~\i - j\ for a XC12 ( |4(a)| | and a YC8 ( J4(b) I 
sample along one lattice axis with exponential and power law 
fits. An x~ 4 line is shown as guide to the eye. 



and depends only on the total quantum dimension D as 
7 = log 2 (D) [731 H3|- In our mappings, DMRG gives 
direct access to density matrices of cylinder slices. We 
calculate S a for cylinders of fixed c and extrapolate in 
L _1 to L — > 00; a linear extrapolation in c — > yields 7. 
Results are ID mapping independent. Wc show interme- 
diate values of a (Fig. [6]), which all show a clearly finite 
value of 7, with a value very consistent with 7 = 1; large- 
a results agree. Small-a results are unreliable, as DMRG 
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FIG. 5. Log-linear plot of the absolute value of the chiral 
correlation function (Co, C x ) = (S io ■ (S jo x S ho ) ■ S i;c ■ (Sj x x 
ofc^)) versus the distance x — |An — A^l along a lattice axis 
for a 196-site YC8 sample with exponential and power law 
fits. 
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FIG. 6. Renyi entropies So, of infinitely long cylinders for 
various a versus circumference c, extrapolated to c = 0. The 
negative intercept is the topological entanglement entropy 7. 

does not capture the tail of the spectrum of p properly, 
but also point to a finite value of 7, hence a topolog- 
ical ground state. The quantum dimension is D = 2, 
excluding chiral spin liquids (7 = 1/2 or D = V? [Z2]). 
Rigorously, DMRG only provides a lower bound on D 
[5D], but the bound is essentially exact as DMRG is a 
method with low entanglement bias [BTj . 

Conclusion. -Through a combination of a large num- 
ber of DMRG states, large samples with small finite 
size effect, and the use of the SU(2) symmetry of the 
kagome model, we have been able to corroborate ear- 
lier evidence for a QSL as opposed to a VBC, due to 
energetic considerations and complete absence of break- 
ing of space group invariance, although DMRG should 
be biased towards VBC due to its low-entanglement na- 
ture and the use of OBC. On the basis of the numerical 
evidence (spin gap, structure factor, spin, dimer and chi- 
ral correlations, topological entanglement entropy) nu- 



merous QSL proposals can be ruled out for the kagome 
system. On the system sizes reached, the spin gap is 
very robust and essentially size- independent, ruling out 
all proposals for gapless spin liquids, consistent with the 
exponential decay of correlators. Individual gapless QSL 
proposals make other predictions not supported by nu- 
merical data, e.g. the static spin structure factor [2"5] , 
Another strong observation is the very rapid decay of 
chiral correlations, ruling out proposals related to chiral 
QSL. The third strong observation is finite topological 
entanglement, which implies a topologically degenerate 
ground state for the kagome system. For quantum di- 
mension 2, as found here, we have in principle, for a 
time-reversal invariant ground state, a choice between 
a Z 2 phase and a double-semion phase [SU [S3]. A Z 2 
QSL emerges straightforwardly in effective field theories 
of the kagome model as a mean-field phase stable un- 
der quantum fluctuations breaking a U(l) gauge sym- 
metry down to Z 2 due to a Higgs mechanism [84], and 
microscopically a resonating valence bond state formed 
from nearest-neighbor Rokhsar-Kivelson dimer coverings 
of the kagome lattice directly leads to a Z 2 QSL [SSJ [SB] 
albeit for a variational energy far from the ground state 
energy. The concentration of weight of the structure fac- 
tor at the hexagonal Brillouin zone edge with shallow 
maxima at the corners would also point to the Z 2 QSL as 
proposed by Sachdev[2"7]. and a Z 2 QSL is also consistent 
with all other numerical findings. All this provides strong 
evidence for the Z 2 QSL, whereas to our knowledge, no 
plausible scenario for the emergence of a double-semion 
phase in the KAFM has been discovered so far, making it 
unplausible, but of course not impossible. An analysis of 
the degenerate ground state manifold as proposed in [SO] , 
not possible with our data, would settle the issue. Even 
if the answer provided final evidence for a Z 2 QSL, many 
questions regarding the detailed microscopic structure of 
the ground state wave functions and the precise nature 
of the Z 2 QSL would remain for future research. 
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Note added. -Recently, we became aware of Ref. [SI] 
which calculates topological entanglement entropy from 
von Neumann entropy for a next-nearest neighbor modifi- 
cation of the KAFM, perfectly consistent with our results 
of D = 2 for the KAFM itself. 
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(a)XC8 system, 38 sites (b) YC6 system, 39 sites 



FIG. 7. (Color online) We map the two-dimensional system 
to a one-dimensional chain using two different mappings, one 
for aligning the lattice to the X-axis (a) and one for aligning 
it to the Y-axis (b) with periodic (open) boundary conditions 
in the vertical (horizontal) direction. The red broad line rep- 
resents the one-dimensional chain. 

Mapping two-dimensional kagome lattices to one- 
dimensional chains for DMRG treatment. As DMRG is 
a one-dimensional method, the two-dimensional kagome 
lattice on cylinders and tori has to be mapped to a one- 
dimensional chain with long-ranged interactions. There 
are multiple (in fact, combinatorially many) ways to map 
cylinders and tori to one-dimensional systems, however, 
ideally they keep interactions as short-ranged and the re- 
sulting path as regular as possible. Out of a large variety 
we tested we choose the two ways to map the kagome 
lattice to chains (Fig. [7| that show the fastest conver- 
gence of energy in DMRG runs and label these either as 
X-cylinders (XC) or Y-cylinders (YC) depending on the 
lattice axes' alignment. In this notation, YC6 denotes a 
cylindrical system where one of the three lattice axes is 
aligned with the y-axis and a circumference c of six lattice 
spacings. For XC systems (alignment of one of the lattice 
axes with the ir-axis) the circumference is measured in- 
stead in units of v3/2 times the lattice spacing, so that 
e.g. the XC12 has a circumference of c = 6^/3 « 10.4 




FIG. 8. (Color online) Path for a torus in DMRG 



(a) Ground state for 5 = 0. 



(a) This snapshot of a not yet converged calculation 
(insufficient number of DMRG states) shows bond 
energy patterns that break the translational invariancc. 




(b) Ground state for 5 = 1. 




(c) Difference in bond energies between the 5 ■ 
the 5 = 1 samples. 



and 



FIG. 9. (Color online) Bond and triangle energies for the 
ground state and lowest triplet excitation of a 74-site XC8 
sample. In panels (a) and (b), the line width is proportional to 
| (Si • Sj) — eo| where eo is the mean bulk energy. Green bonds 
denote {Si ■ Sj) < eo, red (dotted) bonds denote {Si ■ Sj) > e . 
The triangle color (pattern) and intensity correspond to the 
deviation of the sum of the bond energies on the three triangle 
bonds from the mean 3eo, where the green (hatched) triangles 
denote a lower value, i.e. ^ e^ < 3eo. In panel (c), bond 
energies of the lowest S — state are subtracted from those 
of the lowest S = 1 state. The line width is proportional to 
the absolute value of the energy difference, green (hatched) 
lines correspond to positive and red (dotted) lines to negative 
energy differences. 



lattice spacings. In the case of tori, which we consid- 
ered mainly for reference purposes, only a single path 
was retained (Fig. [8]). It is worthwhile to point out the 
path independence of results: where we consider the same 
cylinders as [ST] , results do agree although they used yet 
another mapping. 




(b) This snapshot of a well converged calculation 
(sufficient number of DMRG states and sweeps) shows 
no pattern in the bond energies except for edge effects. 
In the bulk, a spin liquid state without breaking of 
translational invariance emerges. 

FIG. 10. Visualization of the energy per bond for two snap- 
shots in an iterative DMRG ground state calculation. The 
bond line width corresponds to the absolute value of the 
bond energy; the sign is negative (antiferromagnetic) for blue 
bonds, positive (ferromagnetic) for all red bonds, of which 
there are a few towards the edge. 



Identification of bulk vs. boundary excitations. To rule 
out boundary excitations in the lowest S = 1 state, we 
examine the difference in bond energies for the lowest 
lying states in the two spin sectors 5 = and 5=1, 
finding no significant difference at the boundaries but a 
visible change in the bulk (Fig. [9]). 

Supplementary information on ground state properties. 
In order to exclude a valence bond crystal more rigor- 
ously, we consider the bond energies (nearest-neighbor 
correlators) where a valence bond crystal would exhibit 
a frozen pattern of different bond energies. We do not 
observe this for any of our ground states (see Fig. |9ja)). 
Interestingly, it turns out that we can see this frozen 



A 



pattern in unconverged wave functions (Fig 
further increase of the number of kept DMRG states and 
continued sweeping makes these patterns vanish in the 
bulk (Fig. [10(b)] ). The presence of these frozen bond pat- 
terns hence is a distinguishing feature of an insufficiently 
converged wave function as it disappears upon lowering 
the energy and approaching the true ground state, where 
the bond energies only show deviations from the average 
at the cylinder's edges (Fig. [9| . DMRG - similar to other 
tensor network methods such as PEPS and MERA - has 
a low-entanglement bias, because the underlying matrix 
product states structure can only capture entanglement 
up to a strength roughly logarithmic in the number of 
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(a) 



(b) 



(c) 



FIG. 11. (Color online) Snapshots ol the resonance pattern 
for a 150-site YC8 sample, using the same nomenclature as 
Fig. [9] Line widths correspond to the deviation of the bond 
energy from the mean bulk energy eo; triangle color and in- 
tensity show the deviation of the sum of the three triangle 
bond energies from the bulk average 3eo. In each case, in- 
teractions on certain bonds (highlighted by color) have been 
enhanced, in 11(a) a six-site hexag on, in 11(b) an eight-site 
diamond by 0.001 each. In |ll(c)| the interaction strength 
of every second vertical bond was alternatingly changed by 
±0.5%, i.e. every fourth bond was strenghtened. The sur- 
rounding dimers arose in response to these changes, with the 
response increasing from 11(a) to 11(b) to 11(c) 



DMRG states: for an insufficient number of ansatz states, 
DMRG will therefore among states of similar energy pre- 
fer those of low entanglement, in our case valence bond 
crystals compared to quantum spin liquids. 

To elicit additional information on the spin liquid state, 
we strengthen selectively the interaction on various pat- 
terns on some bonds, namely on a hexagon and on a 
diamond pattern and check whether this is taken up by 



the ground state structure (Fig. 11 a) and (b)). In agree- 
ment with the U(l) DMRG calculation of [S3, we find 
in the SU(2) DMRG calculation that strengthening the 
interactions on the diamond pattern elicits the strongest 
response in the bond energies. Agreement is also ob- 
tained for modulating a pattern of every second vertical 



bond (Fig. 11 (c)), which finds an even stronger response; 
this was considered in [57] as evidence that the ground 
state of the kagome model arises from melting a valence 
bond state exhibiting a similar bond pattern. 

As an additional check for preferred orderings, we also 
consider spin-spin correlations in real space (Fig. 12 ) 



where lattice symmetry breaking orderings would show 
up as stronger correlations in certain directions. While 
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FIG. 12. (Color online) Spin-spin correlations in the ground 
state of a 196-site YC8 sample system. The diameter of the 
circles (diamonds) is proportional to the absolute value of 
the spin-spin correlation with the central reference site (black 
square). Blue circles (red diamonds) denote positive (neg- 
ative) correlations. Nearest-neighbor correlations have been 
left our for clarity. The lattice is drawn as a guide for the 
eyes. 



we do not observe any signs for a valence bond crystal 
in the ground state, we see the band-like structure of the 
spin-spin correlations that was reported by Lauchli et 
qL[47] for tori. These pronounced staggered correlations 
along selected loops wrapping around the sample are ar- 
tifacts of the periodic boundary conditions and disappear 
for large circumferences. In agreement with expectations 
for a topological Z2 QSL we also observe the forming of 
band-like structures in the bond energies for cylinders 
with an odd number of sites (not shown). 



